Validating a molecular clock for nudibranchs—No fossils to the rescue

Abstract Time calibrated phylogenies are typically reconstructed with fossil information but for soft‐bodied marine invertebrates that lack hard parts, a fossil record is lacking. In these cases, biogeographic calibrations or the rates of divergence for related taxa are often used. Although nudibranch phylogenies have advanced with the input of molecular data, no study has derived a divergence rate for this diverse group of invertebrates. Here, we use an updated closure date for the Isthmus of Panama (2.8 Ma) to derive the first divergence rates for chromodorid nudibranchs using multigene data from a geminate pair with broad phylogeographic sampling. Examining the species Chromolaichma sedna (Marcus & Marcus, 1967), we uncover deep divergences among eastern Pacific and western Atlantic clades and we erect a new species designation for the latter (Chromolaichma hemera sp. nov.). Next, we discover extensive phylogeographic structure within C. hemera sp. nov. sensu lato, thereby refuting the hypothesis of a recent introduction. Lastly, we derive divergence rates for mitochondrial and nuclear loci that exceed known rates for other gastropods and we highlight significant rate heterogeneity both among markers and taxa. Together, these findings improve understanding of nudibranch systematics and provide rates useful to apply to divergence scenarios in this diverse group.

calibrations, such as the Isthmus of Panama or the closure of the Bering Strait (e.g.Coppard et al., 2013;Hallas et al., 2016;Knowlton & Weigt, 1998;Loeza-Quintana et al., 2018), but these approaches carry a number of assumptions.Most importantly, assigning dates to these events can be difficult to do with precision since dates often span a large period of time (Ho et al., 2015).However, recent work by O'Dea et al. (2016) demonstrated a gradual closure with the Isthmus of Panama forming 2.8 million years ago (Ma), a date that is most appropriate for shallow-water organisms.Despite these advances, biogeographic calibrations can still be problematic when putative geminate species pairs are erroneously assigned and not recovered as sister-species (Marko et al., 2015).
Soft-bodied sea slugs such as nudibranchs, lack diagnostic hard parts large enough (and informative enough) to fossilize and this dearth of fossil data has hindered molecular clock calibration in the group.Although previous studies have calibrated nudibranch phylogenies with few geological data points (Hallas et al., 2016) or rates from other taxonomic groups (Lindsay et al., 2016), a lineage-specific rate is lacking for nudibranchs.Chromolaichma sedna, a chromodorid nudibranch, is known on both sides of the Isthmus of Panama, but was first reported in the western Atlantic only in the 1980s (Bertsch, 1988).A subsequent surge in distribution records for this species in the western Atlantic in the 1990s led to support for a hypothesis of recent introduction in the Caribbean Sea (Goodheart et al., 2016).However, the remarkable conservation of external features in both eastern Pacific and western Atlantic populations (Bertsch, 1988), along with results from a broadly-sampled molecular phylogeny of the Chromodorididae (Johnson & Gosliner, 2012), indicate that these allopatric populations are indeed conspecific and constitute a geminate pair.Investigation of the systematic relationship between these two geographic areas is critical to clarifying their taxonomy.Furthermore, the continued use of few diagnostic loci for nudibranch systematics and phylogenetics (e.g.Jung et al., 2020;Layton et al., 2018;Soong et al., 2020) justifies the need for deriving reliable evolutionary rates for these loci and taxa.
This study seeks to estimate rates of divergence for chromodorid slug species, based on a geminate pair with broad phylogeographic sampling, using the emergence of the Isthmus of Panama as a calibration point.This will provide the first estimate of divergence rates for fossil-less, chromodorid nudibranchs.Next, we address the hypothesis that C. sedna is a recent introduction to the western Atlantic by investigating phylogeographic structure in the region.Finally, we revise the taxonomic status of the western Atlantic Chromolaichma 'sedna'.

| Specimen collection and DNA sequencing
A total of 20 specimens of C. 'sedna' were collected by hand on SCUBA at depths of 5-15 m at multiple sites in the Caribbean Sea and eastern Pacific Ocean, including near the type locality in the Sea of Cortez (Figure 1, Table 1).Two additional specimens were provided by Florida Natural History Museum and Museum of Comparative Zoology and sequence data for three individuals were sourced from Genbank (Table 1).All specimens were preserved in 96% ethanol in preparation for molecular analyses.Genomic DNA was extracted from a sample of the foot tissue using a DNeasy blood & tissue kit (Qiagen) following manufacturer's guidelines.PCR amplicons were purified using either Exo-SapIT or gel purification.Two mitochondrial genes, cytochrome c oxidase subunit I (COI) and 16S ribosomal RNA (16S), and one nuclear gene, adenine nucleotide translocase (ANT), were targeted for sequencing.For PCR and sequencing, Folmer et al. (1994) primers were used for COI, Palumbi et al. (1991) primers were used for 16S, and Audzijonyte and Vrijenhoek (2010) primers were used for ANT.Amplicons were bidirectionally sequenced on an Applied Biosystems 3730 capillary sequencer.Forward and reverse sequences were edited and assembled in Sequencher (Gene Codes Corporation) and Geneious Prime 2021 (http:// www.genei ous.com/ ) and aligned with MAFFT (Katoh et al., 2002).A second alignment was trimmed to equal length in Geneious Prime 2021 for haplotype network construction.

| Phylogenetic analyses and species delimitation
Both individual gene and concatenated gene datasets were used for constructing maximum likelihood (ML) trees in RAxML (Stamatakis, 2006) implemented in the raxmlGUI v2.0 (Edler et al., 2020).Trees were built with a GTR + G model and 1000 bootstrap replicates with Doriprismatica atromarginata as an outgroup for the concatenated dataset, COI and 16S, and with Tyrinna evelinae as an outgroup for ANT.Total branch length (TBL) of each gene tree was summed and used to derive a mutation rate ratio for mitochondrial (μm) and nuclear loci (μn) following Duda (2021).For μm, an average TBL of both genes was used.Both the Assemble Species by Automatic Partition (ASAP) algorithm (Puillandre et al., 2021) and the Bayesian bPTPmodel (Zhang et al., 2013) were used to partition the COI dataset into unique genetic clusters using default parameters.A minimum spanning haplotype network was generated with trimmed COI data using the haploNet function in the pegas package in R (Paradis, 2010).

| Estimating divergence rates
Pairwise distances were calculated for COI, 16S, and ANT in MEGAX (Kumar et al., 2018) using a maximum composite likelihood model and pairwise deletion.To calculate a rate of divergence for this species, we used the minimum sequence divergence between geminate pairs and a closing date of the Isthmus of Panama of 2.8 mya (O'Dea et al., 2016).We chose to use the more conservative minimum sequence divergence rather than mean divergence to prevent the inclusion of accumulated intraspecific changes.Rates of divergence were calculated for all genes.The substitution rate is presented as half of the divergence rate.To determine whether the data were evolving in a clock-like fashion, a maximum likelihood clock test was conducted in MEGAX (Kumar et al., 2018) using a Tamura-Nei model and rate heterogeneity for complete COI, 16S and ANT sequences.Divergence rates from other geminate pairs of marine invertebrates were retrieved from the literature for comparative purposes, along with relevant life history information.Life history information, including lifespan, reproductive mode and larval dispersal, were obtained for the genus where possible but may represent broader taxonomic patterns in some instances.

| Deep divergence within and among ocean basins
A total 22 COI, 22 16S and 19 ANT sequences were generated for C. 'sedna' in this study, with 25 total sequences of COI and 16S including those sourced from GenBank.Data for these specimens and GenBank accession numbers are provided in Table 1.The final COI alignment was 658 bp in length, the 16S alignment was 459 bp and the ANT alignment was 553 bp.A trimmed COI alignment of 590 bp was used for the haplotype network.Mean intraspecific divergence of C. 'sedna' was 8.7% (range 0.0%-15.7%)at COI, 3.9% (range 0.0%-8.1%)at 16S and 0.3% (range 0.0%-0.7%)at ANT (Tables S1-S3).
The COI + 16S + ANT ML phylogeny showed that western Atlantic and eastern Pacific C. 'sedna' formed separate clades, with more structure in the western Atlantic than in the eastern Pacific (Figure 2a).These same patterns were recovered in individual gene trees (Figure S1).Three western Atlantic clades were strongly supported in the combined phylogeny (BP > 90), although support at interior nodes was weak (Figure 2a).These results were reinforced by the species delimitation analyses that partitioned western Atlantic individuals into three clusters, corresponding to geographic region (Bahamas, Panama, and Florida) (Figure 2a).The minimum divergence between individuals in the western Atlantic and eastern Pacific was 13.6% at COI, 5.7% at 16S and 0.5% at ANT.Within the western Atlantic clade mean intraspecific divergence at COI was 4.0%, with a mean divergence of 5.4% between Florida and Bahamas, 6.5% between Bahamas and Panama, and 6.5% between Florida and Panama.The minimum sequence divergence at COI among these clades were 5.4%, 6.3%, and 6.2%, respectively.Within the eastern Pacific clade mean intraspecific divergence at COI was 0.45%.

| Rate variation
Given a minimum sequence divergence of 13.6% among C. 'sedna' individuals on either side of the Isthmus of Panama, and an estimated age of separation of 2.8 million years, COI is expected to diverge at a rate of 4.9% per million years.For 16S, with a minimum sequence divergence of 5.7%, this rate was 2.0% per million years.For ANT, with a minimum sequence divergence of 0.5%, this rate was 0.2% per million years.Divergence rates of 4.9%, 2.0%, and 0.2% correspond to substitution rates of 2.5% for COI, 1.0% for 16S, and 0.1% for ANT.
COI and 16S rates are 25 and 10 times higher than the ANT rates, respectively.This variation in divergence rate also underpins differences in pairwise distance among loci (Figure S2).Total branch length of the COI, 16S, and ANT trees was 0.299, 0.136, and 0.006, respectively, with a μm/μn value of 39.5.The clock test was insignificant for COI, 16S, and ANT (p = .9,p = .7,p = .06),indicating that the null hypothesis of equal evolutionary rates throughout the tree could not be rejected (Table 2).
The central rachidian tooth is a small, reduced, triangular tooth.The innermost tooth in each row had 2-4 denticles on the inner side and 3-5 outer side denticles.The rest of the teeth in each half row have 4-8 denticles on the outer side, but these became smooth for the outermost 10-20 teeth.The jaw rodlets were mostly bifid.
The reproductive system of these two specimens were examined by Bertsch (1988), who stated that the morphology of the penis, vas deferens, and the arrangement and relative size of the vagina, insemination duct, bursa copulatrix and the receptaculum seminis were all identical to the holotype of C. sedna (Marcus & Marcus, 1967).

| Remarks
Chromolaichma was only recently resurrected by Matsuda and Gosliner (2018)  Morphological conservation between C. hemera sp.nov.and C. sedna was noted by Bertsch (1988) who described them as nearly identical, noting slight differences in background color (not as white as C. sedna) and with possibly more crenulations in the mantle (Figure 4).Chromolaichma hemera sp.nov.also exhibits a wider yellow submarginal band than C. sedna (Figure 4).Bertsch noted that the reproductive system was identical to that illustrated in the original description of C. sedna (Marcus & Marcus, 1967) and that the radular morphology matched the holotype, and fell within the range for specimens described from the Gulf of California (Bertsch, 1978).Overall, he did not consider these slight morpho-   Do et al., 2020;Epstein et al., 2019;Matsuda & Gosliner, 2018), but the nudibranch fauna of the western Atlantic has been less well studied than the Indo-Pacific.
Our study also uncovered significant population structure within western Atlantic C. hemera sp.nov., as evidenced by deep haplotypic divergence among populations, negating the hypothesis that this species was recently introduced to the area through ship transport (e.g.Goodheart et al., 2016;Humann & Deloach, 2002).Bertsch (1988) mistakenly reported the location for a specimen shared with him by Lyons when the first tropical western Atlantic specimens were reported (see Lyons, 1989 for correction).Thus, Bertsch incorrectly reported that all known specimens were from a restricted locality, initiating the idea that they may represent an introduction, although not explicitly stating that.The earliest known specimens in the western Atlantic were identified from the Bahamas in 1978 (Lyons, 1989) nov.sensu lato).Conversely, we uncovered a lack of structure within the eastern Pacific clade.Although surprising, this has also been shown in Tripneustes sea urchins (Lessios et al., 2003) and the pronghorn spiny lobster (Panulirus penicillatus) (Chow et al., 2011) and it may be due to differences in oceanographic currents, habitat complexity and glacial sea level changes in the region.
There are several possible explanations for the deep divergence among western Atlantic clades.First, these deep divergences might reflect missing taxa given the high extinction rate in the region following the closure of the Isthmus (Lessios, 2008;Marko et al., 2015), or they might reflect rapid diversification following the split (Thacker, 2017)

| Rate heterogeneity among geminate lineages
Rate heterogeneity has been invoked to explain variation in sequence divergence observed across marine geminate pairs (Marko, 2002), but this is the first study to summarize these rates.Here, we found significant rate heterogeneity among marine invertebrate genera that share a similar lifespan, reproductive mode and dispersal potential, suggesting that life history traits do not drive differences in mitochondrial rates.These findings align with recent work by Saclier et al. (2018), who found that while nuclear substitution rates were linked to generation time in isopods, mitochondrial substitution rates were independent of life history traits.Some of the seminal work linking mutation rates to life history traits was focused on only a few taxa and loci (e.g.Müller & Albach, 2010), thereby limiting our understanding of rate variation across the genomes of most species.In fact, we uncovered significant differences in nuclear and mitochondrial substitution and mutation rates, as evidenced by an elevated μm/μn ratio, similar to findings in other gastropods (Duda, 2021).
Given the evidence for rate variation among both taxa and loci, future work will look to reassess substitution rates using transcriptome data from recent phylogenomic studies (Layton et al., 2020;Moles & Giribet, 2021).However, because COI, 16S, and ANT remain some of the most commonly used genes in species-level nudibranch investigations, the lineage-specific rates presented here will help improve understanding of divergence scenarios in this diverse group.

AUTH
The same patterns were also recovered in the haplotype network analysis, with different ocean basins showing significant F I G U R E 1 Sampling locations for 25 Chromolaichma specimens collected from either side of the Isthmus of Panama.Sample numbers appear next to locations.Images derived from iNaturalist (on left, photographer: S. & J. Johnson) and Museum of Comparative Zoology (on right, photographer: G. Giribet).TA B L E 1 Specimen information for 25 Chromolaichma 'sedna' samples collected in the Caribbean Sea and Pacific Ocean for this study, along with outgroup taxa.GenBank accession numbers are provided for COI, 16S and ANT.An asterisk marks samples that derive from GenBank and a hyphen denotes missing information.differentiation, as well as regions within the western Atlantic showing deep haplotype divergence and no shared haplotypes (Figure2b).In total, 8 and 9 unique haplotypes were identified from 12 and 13 western Atlantic and eastern Pacific individuals, respectively (Figure2b).Within the Pacific, all haplotypes were closely related to the central haplotype in Mexico and the Sea of Cortez.

3. 3
.5 | Description Body soft, elongate and translucent white, with some specimens (including the holotype) with some orange-red dusting on the dorsum and sides of foot.Mantle edge usually crenulated.Tail extends beyond mantle.Mantle with yellow marginal band, and thinner red submarginal band.Rarely, an opaque band sits inside the red submarginal band.Similar bands of color appear on the foot margin, although the red submarginal band can be very thin or absent.Rhinophore stalks, translucent white.Lower part of the rhinophore clubs are translucent white, with anterior top third lamellae red.For most specimens, posteriorly, the stem that the rhinophore lamellae are attached to is red, and anteriorly, white.Lamellae evenly red around the upper part of rhinophores.Some specimens show lateral rhinophoral lamellae as white.Simple pinnate gills form a double spiral, and wriggle rhythmically.Large specimens may be bipinnate at tips.Distal tips of outer rachis and lamellae are red.
logical differences to warrant a new species, but the addition of molecular data that demonstrates deep genetic divergences does support the establishment of a new species for the western Atlantic and Caribbean specimens.Molecular data also demonstrates deep divergences among western Atlantic and Caribbean clades, which may represent additional cryptic species.A lack of geographically intermediate samples and additional morphological data prevents their description here, so for now, we refer to all lineages in the western Atlantic and Caribbean clades as C. hemera sp.nov.sensu lato.Upon interrogation of publicly available images, there are notable differences in rhinophore colouration between individuals in the lesser Antilles and the rest of the Caribbean, but additional work is needed to verify this preliminary finding.Finally, the name C. sedna is now restricted to specimens in the Gulf of California and the eastern Pacific.
, only 11 years after the description of the eastern Pacific species C. sedna.The deep molecular divergences among tropical western Atlantic populations do not support a recent introduction scenario and instead may represent a cryptic species complex (C.hemera sp. O R CO NTR I B UTI O N S Kara K. S. Layton: Conceptualization (supporting); formal analysis (lead); investigation (equal); validation (equal); visualization (equal); writing -original draft (equal); writing -review and editing (equal).Nerida G. Wilson: Conceptualization (lead); formal analysis (supporting); funding acquisition (lead); investigation (equal); resources (lead); validation (equal); visualization (equal); writing -original draft (equal); writing -review and editing (equal).ACK N OWLED G M ENTS Greg Rouse (UCSD-SIO), Pat Krug (Cal State LA), Liz Borda (Texas A&M), and Laeticia Serrano assisted with specimen collection and provided samples.John Slapscinsky and Mandy Bemis facilitated specimen loans from the Florida Museum of Natural History, Charlotte Seid facilitated loans from Scripps Institution of Oceanography, Benthic Invertebrates Collection and Jennifer Trimble and Emily Blank facilitated a specimen loan from the Museum of Comparative Zoology (Harvard).New collections for this study were made in the Gulf of California, Mexico under permit DGOPA.06754.240609.2003 to Carlos Ortiz, Universidad Autónoma de Baja Californias Sur.New collections in Florida The number of species per taxonomic group was sourced from WoRMS, based on accepted records only, and are presented for both the genus and family (in brackets) where applicable.Information about th lifespan, reproductive mode, and larval dispersal were obtained for the genus where possible but may represent broader taxonomic patterns in some instances.Arrows indicate high or low larval dispersal potential and "y" represents the year.